##  Randomization inference for EXPERIMENT 2 for sponsorship paper
rm(list = ls())

library(acs)
library(RCurl)
library(RJSONIO)


setwd("~/Dropbox/Current projects/kern-pietryka-crabtree/data/")
dat <- read.csv("sponsorship_2_January+10%2C+2018_12.58.csv", header = TRUE)
dat <- dat[-c(1:2),]

dim(dat)
colnames(dat)


##  consistency check: only real responses are included
table(dat$Status)


##  drop responses that were not completed
sel <- dat$Finished == "True"
table(sel)
dat <- dat[sel,]


##  drop respondents who don't give consent
sel <- dat$Q2 == "Yes"
table(sel) 
dat <- dat[sel,]


##  any respondents who answer multiple times?
stopifnot(length(dat$mTurkCode) == length(unique(dat$mTurkCode)))


##  drop respondents based on lat/long
latlong2fips <- function(latitude, longitude) {
  url <- "http://data.fcc.gov/api/block/find?format=json&latitude=%f&longitude=%f"
  url <- sprintf(url, latitude, longitude)
  json <- RCurl::getURL(url)
  json <- RJSONIO::fromJSON(json)
  as.character(json$County['FIPS'])
}
key <- api.key.install("f66712573be9fdfd1264a2d99bc0356a4b3edaf7")



###  this needs to be fixed to work properly; new web address

latlong2fips <- 
function(latitude, longitude) {
  url <- "https://geo.fcc.gov/api/census/block/find?latitude=%f&longitude=%f&format=json"
  url <- sprintf(url, latitude, longitude)
  json <- RCurl::getURL(url)
  json <- RJSONIO::fromJSON(json)
  as.character(json$County['FIPS'])
}



dat$state <- NA

m <- 1
for(m in 1:nrow(dat))	{

	temp <- latlong2fips(latitude = as.numeric(as.character(dat$LocationLatitude[m])),
						 longitude = as.numeric(as.character(dat$LocationLongitude[m])))

if(is.null(temp))	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(is.na(temp)) 	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(temp == "NULL")	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }

	sel <- which(fips.state[1:51,1] == as.numeric(substr(temp,1,2)))
	temp <- fips.state[sel,2]

	if(length(temp) == 0)	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
								as.numeric(as.character(dat$LocationLongitude[m])),"\n");
								next() }
	dat$state[m] <- temp

	if(m %% 100 == 0) cat(m, "\n")

						}

##  drop respondents not from the 50 states
sel <- is.na(dat$state)
table(sel)
dat <- dat[!sel,]
dim(dat)

table(dat$state)


##  Deal with identical lat/long coordinates
##  Drop respondent if previous respondent with identical lat/long exists in data
dat$latlong <- paste(dat$LocationLatitude, dat$LocationLongitude, sep = "/")

dat$end.date <- strptime(dat[,2], format = "%Y-%m-%d %H:%M:%S")
dat <- dat[order(dat$end.date),]

##  consistency check
stopifnot(all(dat$end.date == sort(dat$end.date)))

dat$repeated <- NA
dat$repeated[1] <- FALSE


m <- 2
for(m in 2:nrow(dat))	{

	if(!any(dat$latlong[m] == dat$latlong[1:(m-1)])) { dat$repeated[m] <- FALSE } else
	{ dat$repeated[m] <- TRUE }

						}
table(dat$repeated)
sum(is.na(dat$repeated))
dat <- dat[!dat$repeated,]



##  baseline questions
dat$show <- dat$Q37
dat$show <- droplevels(dat$show)
table(dat$show)
dat$show.n <- dat$show
levels(dat$show.n) <- c(3,5,4,1,2)
table(dat$show,dat$show.n)


dat$actor <- dat$Q38
dat$actor <- droplevels(dat$actor)
table(dat$actor)
dat$actor.n <- dat$actor
levels(dat$actor.n) <- c(5,3,4,2,1)
table(dat$actor,dat$actor.n)


dat$attention <- dat$Q39
dat$attention <- droplevels(dat$attention)
table(dat$attention)
dat$attention.n <- dat$attention
levels(dat$attention.n) <- c(3,5,4,1,2)
table(dat$attention,dat$attention.n)


dat$index <-	as.numeric(as.character(dat$show.n)) + 
				as.numeric(as.character(dat$actor.n)) + 
				as.numeric(as.character(dat$attention.n))

barplot(table(dat$index),
col = c(rep("lightgreen",4),rep("grey",2),rep("cornflowerblue",7)) ,
main = "Self-monitoring scale", xlab = "", cex.lab = 1.2)



dat$Tr <- NA
dat$liberals <- NA
dat$conservatives <- NA
dat$muslims <- NA
dat$welfare <- NA
dat$atheists <- NA
dat$feminists <- NA


m <- 1
for(m in 1:nrow(dat))	{

	temp <- 			c(	as.character(dat$Q42_1[m]),
							as.character(dat$Q61_1[m]),
							as.character(dat$Q48_1[m]))

	dat$liberals[m] <- temp[temp != ""]
	dat$Tr[m] <- which(temp != "")	##  order: no logo, UC Berkeley, Liberty U

	temp <- 			c(	as.character(dat$Q43_1[m]),
							as.character(dat$Q62_1[m]),
							as.character(dat$Q49_1[m]))
	dat$conservatives[m] <- temp[temp != ""]


	temp <- 			c(	as.character(dat$Q44_1[m]),
							as.character(dat$Q63_1[m]),
							as.character(dat$Q50_1[m]))
	dat$muslims[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q45_1[m]),
							as.character(dat$Q64_1[m]),
							as.character(dat$Q51_1[m]))
	dat$welfare[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q46_1[m]),
							as.character(dat$Q65_1[m]),
							as.character(dat$Q52_1[m]))
	dat$atheists[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q40_1[m]),
							as.character(dat$Q66_1[m]),
							as.character(dat$Q53_1[m]))
	dat$feminists[m] <- temp[temp != ""]

						}

table(dat$Tr)

##  test whether random assignment is as good as random
chisq.test(table(dat$Tr), simulate.p.value = TRUE, B = 10000)


dat$liberals <- as.numeric(as.character(dat$liberals))
dat$conservatives <- as.numeric(as.character(dat$conservatives))
dat$muslims <- as.numeric(as.character(dat$muslims))
dat$welfare <- as.numeric(as.character(dat$welfare))
dat$atheists <- as.numeric(as.character(dat$atheists))
dat$feminists <- as.numeric(as.character(dat$feminists))

##  we need to reverse code conservative
dat$conservatives <- abs(dat$conservatives-100)

# save.image("Experiment 2 intermediate results.RData")
# load("Experiment 2 intermediate results.RData")



##  randomization inference
s <- 1000000

##  compute sample test statistics for all six outcomes and all three-way comparisons
##  note that tests are one-sided here
dat$liberals.r <- rank(dat$liberals)
dat$conservatives.r <- rank(dat$conservatives)
dat$muslims.r <- rank(dat$muslims)
dat$welfare.r <- rank(dat$welfare)
dat$atheists.r <- rank(dat$atheists)
dat$feminists.r <- rank(dat$feminists)


##  Order: no logo, UC Berkeley, Liberty U
##  We generate test statistics in such a way that larger values mean rejection of the null
t.stat.1.1v2 <- mean(dat$liberals.r[dat$Tr == 2]) - mean(dat$liberals.r[dat$Tr == 1])
t.stat.1.1v3 <- mean(dat$liberals.r[dat$Tr == 1]) - mean(dat$liberals.r[dat$Tr == 3])
t.stat.1.2v3 <- mean(dat$liberals.r[dat$Tr == 2]) - mean(dat$liberals.r[dat$Tr == 3])

t.stat.2.1v2 <- mean(dat$conservatives.r[dat$Tr == 2]) - mean(dat$conservatives.r[dat$Tr == 1])
t.stat.2.1v3 <- mean(dat$conservatives.r[dat$Tr == 1]) - mean(dat$conservatives.r[dat$Tr == 3])
t.stat.2.2v3 <- mean(dat$conservatives.r[dat$Tr == 2]) - mean(dat$conservatives.r[dat$Tr == 3])

t.stat.3.1v2 <- mean(dat$muslims.r[dat$Tr == 2]) - mean(dat$muslims.r[dat$Tr == 1])
t.stat.3.1v3 <- mean(dat$muslims.r[dat$Tr == 1]) - mean(dat$muslims.r[dat$Tr == 3])
t.stat.3.2v3 <- mean(dat$muslims.r[dat$Tr == 2]) - mean(dat$muslims.r[dat$Tr == 3])

t.stat.4.1v2 <- mean(dat$welfare.r[dat$Tr == 2]) - mean(dat$welfare.r[dat$Tr == 1])
t.stat.4.1v3 <- mean(dat$welfare.r[dat$Tr == 1]) - mean(dat$welfare.r[dat$Tr == 3])
t.stat.4.2v3 <- mean(dat$welfare.r[dat$Tr == 2]) - mean(dat$welfare.r[dat$Tr == 3])

t.stat.5.1v2 <- mean(dat$atheists.r[dat$Tr == 2]) - mean(dat$atheists.r[dat$Tr == 1])
t.stat.5.1v3 <- mean(dat$atheists.r[dat$Tr == 1]) - mean(dat$atheists.r[dat$Tr == 3])
t.stat.5.2v3 <- mean(dat$atheists.r[dat$Tr == 2]) - mean(dat$atheists.r[dat$Tr == 3])

t.stat.6.1v2 <- mean(dat$feminists.r[dat$Tr == 2]) - mean(dat$feminists.r[dat$Tr == 1])
t.stat.6.1v3 <- mean(dat$feminists.r[dat$Tr == 1]) - mean(dat$feminists.r[dat$Tr == 3])
t.stat.6.2v3 <- mean(dat$feminists.r[dat$Tr == 2]) - mean(dat$feminists.r[dat$Tr == 3])


t.stat.rand <- matrix(NA,s,18)


set.seed(1)
for(m in 1:s)	{

	dat$Tr.r <- sample(c(1,2,3), nrow(dat), replace = TRUE)

	t.stat.rand[m,1] <- mean(dat$liberals.r[dat$Tr.r == 2]) - mean(dat$liberals.r[dat$Tr.r == 1])
	t.stat.rand[m,2] <- mean(dat$liberals.r[dat$Tr.r == 1]) - mean(dat$liberals.r[dat$Tr.r == 3])
	t.stat.rand[m,3] <- mean(dat$liberals.r[dat$Tr.r == 2]) - mean(dat$liberals.r[dat$Tr.r == 3])

	t.stat.rand[m,4] <- mean(dat$conservatives.r[dat$Tr.r == 2]) - mean(dat$conservatives.r[dat$Tr.r == 1])
	t.stat.rand[m,5] <- mean(dat$conservatives.r[dat$Tr.r == 1]) - mean(dat$conservatives.r[dat$Tr.r == 3])
	t.stat.rand[m,6] <- mean(dat$conservatives.r[dat$Tr.r == 2]) - mean(dat$conservatives.r[dat$Tr.r == 3])

	t.stat.rand[m,7] <- mean(dat$muslims.r[dat$Tr.r == 2]) - mean(dat$muslims.r[dat$Tr.r == 1])
	t.stat.rand[m,8] <- mean(dat$muslims.r[dat$Tr.r == 1]) - mean(dat$muslims.r[dat$Tr.r == 3])
	t.stat.rand[m,9] <- mean(dat$muslims.r[dat$Tr.r == 2]) - mean(dat$muslims.r[dat$Tr.r == 3])

	t.stat.rand[m,10] <- mean(dat$welfare.r[dat$Tr.r == 2]) - mean(dat$welfare.r[dat$Tr.r == 1])
	t.stat.rand[m,11] <- mean(dat$welfare.r[dat$Tr.r == 1]) - mean(dat$welfare.r[dat$Tr.r == 3])
	t.stat.rand[m,12] <- mean(dat$welfare.r[dat$Tr.r == 2]) - mean(dat$welfare.r[dat$Tr.r == 3])

	t.stat.rand[m,13] <- mean(dat$atheists.r[dat$Tr.r == 2]) - mean(dat$atheists.r[dat$Tr.r == 1])
	t.stat.rand[m,14] <- mean(dat$atheists.r[dat$Tr.r == 1]) - mean(dat$atheists.r[dat$Tr.r == 3])
	t.stat.rand[m,15] <- mean(dat$atheists.r[dat$Tr.r == 2]) - mean(dat$atheists.r[dat$Tr.r == 3])

	t.stat.rand[m,16] <- mean(dat$feminists.r[dat$Tr.r == 2]) - mean(dat$feminists.r[dat$Tr.r == 1])
	t.stat.rand[m,17] <- mean(dat$feminists.r[dat$Tr.r == 1]) - mean(dat$feminists.r[dat$Tr.r == 3])
	t.stat.rand[m,18] <- mean(dat$feminists.r[dat$Tr.r == 2]) - mean(dat$feminists.r[dat$Tr.r == 3])


	if(m %% 1000 == 0)	{ cat(m, "\n") }

				}

# save.image("Experiment 2 results.RData")
# load("Experiment 2 results.RData")



##  p-values for each outcome and each two-way comparison
##  in each case the expectation is that the test statistic is positive
p <- matrix(NA,6,4)
rownames(p) <- c("Liberals","Conservatives","Muslims","People on welfare","Atheists","Feminists")
colnames(p) <- c("No logo vs. UC Berkeley","No logo vs. Liberty U","UC Berkeley vs. Liberty U","Joint")


p[1,1:3] <- 
c(
if(t.stat.1.1v2 > 0) round(mean(t.stat.rand[,1] > t.stat.1.1v2),2) else "*",
if(t.stat.1.1v3 > 0) round(mean(t.stat.rand[,2] > t.stat.1.1v3),2) else "*", 
if(t.stat.1.2v3 > 0) round(mean(t.stat.rand[,3] > t.stat.1.2v3),2) else "*")

p[2,1:3] <- 
c(
if(t.stat.2.1v2 > 0) round(mean(t.stat.rand[,4] > t.stat.2.1v2),2) else "*",
if(t.stat.2.1v3 > 0) round(mean(t.stat.rand[,5] > t.stat.2.1v3),2) else "*",
if(t.stat.2.2v3 > 0) round(mean(t.stat.rand[,6] > t.stat.2.2v3),2) else "*")

p[3,1:3] <- 
c(
if(t.stat.3.1v2 > 0) round(mean(t.stat.rand[,7] > t.stat.3.1v2),2) else "*",
if(t.stat.3.1v3 > 0) round(mean(t.stat.rand[,8] > t.stat.3.1v3),2) else "*",
if(t.stat.3.2v3 > 0) round(mean(t.stat.rand[,9] > t.stat.3.2v3),2) else "*")

p[4,1:3] <- 
c(
if(t.stat.4.1v2 > 0) round(mean(t.stat.rand[,10] > t.stat.4.1v2),2) else "*",
if(t.stat.4.1v3 > 0) round(mean(t.stat.rand[,11] > t.stat.4.1v3),2) else "*",
if(t.stat.4.2v3 > 0) round(mean(t.stat.rand[,12] > t.stat.4.2v3),2) else "*")

p[5,1:3] <- 
c(
if(t.stat.5.1v2 > 0) round(mean(t.stat.rand[,13] > t.stat.5.1v2),2) else "*",
if(t.stat.5.1v3 > 0) round(mean(t.stat.rand[,14] > t.stat.5.1v3),2) else "*",
if(t.stat.5.2v3 > 0) round(mean(t.stat.rand[,15] > t.stat.5.2v3),2) else "*")

p[6,1:3] <- 
c(
if(t.stat.6.1v2 > 0) round(mean(t.stat.rand[,16] > t.stat.6.1v2),2) else "*",
if(t.stat.6.1v3 > 0) round(mean(t.stat.rand[,17] > t.stat.6.1v3),2) else "*",
if(t.stat.6.2v3 > 0) round(mean(t.stat.rand[,18] > t.stat.6.2v3),2) else "*")


##  p-values for family-wise comparisons (within outcomes)
##  We do not use Young's approach because the test is one-sided
p[,4] <- 
round(
c(
mean(apply(t.stat.rand[,1:3],1,max) 	> max(c(t.stat.1.1v2,t.stat.1.1v3,t.stat.1.2v3))),
mean(apply(t.stat.rand[,4:6],1,max) 	> max(c(t.stat.2.1v2,t.stat.2.1v3,t.stat.2.2v3))),
mean(apply(t.stat.rand[,7:9],1,max) 	> max(c(t.stat.3.1v2,t.stat.3.1v3,t.stat.3.2v3))),
mean(apply(t.stat.rand[,10:12],1,max) 	> max(c(t.stat.4.1v2,t.stat.4.1v3,t.stat.4.2v3))),
mean(apply(t.stat.rand[,13:15],1,max) 	> max(c(t.stat.5.1v2,t.stat.5.1v3,t.stat.5.2v3))),
mean(apply(t.stat.rand[,16:18],1,max) 	> max(c(t.stat.6.1v2,t.stat.6.1v3,t.stat.6.2v3)))),2)


##  omnibus p-value
round(
mean(apply(t.stat.rand[,1:18],1,max) > max(c(	t.stat.1.1v2,t.stat.1.1v3,t.stat.1.2v3,
												t.stat.2.1v2,t.stat.2.1v3,t.stat.2.2v3,
												t.stat.3.1v2,t.stat.3.1v3,t.stat.3.2v3,
												t.stat.4.1v2,t.stat.4.1v3,t.stat.4.2v3,
												t.stat.5.1v2,t.stat.5.1v3,t.stat.5.2v3,
												t.stat.6.1v2,t.stat.6.1v3,t.stat.6.2v3))
												),2)

library(Hmisc)
latex(p, file = "")



##
##  Treatment effect heterogeneity
##
quantile(dat$index, probs = c(0,1/3,2/3,1))
dat$low <- dat$index < 7
dat$high <- dat$index > 8

table(dat$low)
table(dat$high)
sum(dat$low == FALSE & dat$high == FALSE)

# order: no logo, UC Berkeley, Liberty U


tau.high <-	mean(dat$liberals.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$liberals.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.1.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.1.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$liberals.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$liberals.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.1.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.1.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$liberals.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$liberals.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$liberals.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.1.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.1.2v3.h <- -(tau.high-tau.low)



tau.high <-	mean(dat$conservatives.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$conservatives.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.2.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.2.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$conservatives.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$conservatives.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.2.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.2.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$conservatives.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$conservatives.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$conservatives.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.2.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.2.2v3.h <- -(tau.high-tau.low)



tau.high <-	mean(dat$muslims.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$muslims.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.3.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.3.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$muslims.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$muslims.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.3.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.3.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$muslims.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$muslims.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$muslims.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.3.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.3.2v3.h <- -(tau.high-tau.low)



tau.high <-	mean(dat$welfare.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$welfare.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.4.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.4.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$welfare.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$welfare.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.4.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.4.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$welfare.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$welfare.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$welfare.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.4.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.4.2v3.h <- -(tau.high-tau.low)



tau.high <-	mean(dat$atheists.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$atheists.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.5.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.5.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$atheists.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$atheists.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.5.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.5.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$atheists.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$atheists.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$atheists.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.5.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.5.2v3.h <- -(tau.high-tau.low)



tau.high <-	mean(dat$feminists.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 1 & dat$high == TRUE])
tau.low <-	mean(dat$feminists.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 1 & dat$low == TRUE])
t.stat.6.1v2.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.6.1v2.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$feminists.r[dat$Tr == 1 & dat$high == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$feminists.r[dat$Tr == 1 & dat$low == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.6.1v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.6.1v3.h <- -(tau.high-tau.low)


tau.high <-	mean(dat$feminists.r[dat$Tr == 2 & dat$high == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 3 & dat$high == TRUE])
tau.low <-	mean(dat$feminists.r[dat$Tr == 2 & dat$low == TRUE]) -
			mean(dat$feminists.r[dat$Tr == 3 & dat$low == TRUE])
t.stat.6.2v3.h <- tau.high-tau.low
if(tau.high > tau.low & sign(tau.high) == -1) t.stat.6.2v3.h <- -(tau.high-tau.low)



s <- 1000000
t.stat.rand.het <- matrix(NA,s,18)


set.seed(1)
m <- 1
for(m in 1:s)	{

	dat$Tr.r <- sample(c(1,2,3), nrow(dat), replace = TRUE)


	tau.high <-	mean(dat$liberals.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$liberals.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,1] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,1] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$liberals.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$liberals.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,2] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,2] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$liberals.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$liberals.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$liberals.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,3] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,3] <- -(tau.high-tau.low)



	tau.high <-	mean(dat$conservatives.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$conservatives.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,4] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,4] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$conservatives.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$conservatives.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,5] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,5] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$conservatives.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$conservatives.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$conservatives.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,6] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,6] <- -(tau.high-tau.low)



	tau.high <-	mean(dat$muslims.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$muslims.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,7] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,7] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$muslims.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$muslims.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,8] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,8] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$muslims.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$muslims.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$muslims.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,9] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,9] <- -(tau.high-tau.low)



	tau.high <-	mean(dat$welfare.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$welfare.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,10] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,10] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$welfare.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$welfare.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,11] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,11] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$welfare.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$welfare.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$welfare.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,12] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,12] <- -(tau.high-tau.low)



	tau.high <-	mean(dat$atheists.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$atheists.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,13] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,13] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$atheists.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$atheists.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,14] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,14] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$atheists.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$atheists.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$atheists.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,15] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,15] <- -(tau.high-tau.low)



	tau.high <-	mean(dat$feminists.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 1 & dat$high == TRUE])
	tau.low <-	mean(dat$feminists.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 1 & dat$low == TRUE])
	t.stat.rand.het[m,16] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,16] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$feminists.r[dat$Tr.r == 1 & dat$high == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$feminists.r[dat$Tr.r == 1 & dat$low == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,17] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,17] <- -(tau.high-tau.low)


	tau.high <-	mean(dat$feminists.r[dat$Tr.r == 2 & dat$high == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 3 & dat$high == TRUE])
	tau.low <-	mean(dat$feminists.r[dat$Tr.r == 2 & dat$low == TRUE]) -
				mean(dat$feminists.r[dat$Tr.r == 3 & dat$low == TRUE])
	t.stat.rand.het[m,18] <- tau.high-tau.low
	if(tau.high > tau.low & sign(tau.high) == -1) t.stat.rand.het[m,18] <- -(tau.high-tau.low)


	if(m %% 1000 == 0)	{ cat(m, "\n") }

				}

# save.image("Experiment 2 results het.RData")
# load("Experiment 2 results het.RData")



##  p-values for each outcome and each two-way comparison
##  in each case the expectation is that the test statistic is positive
p.het <- matrix(NA,6,4)
rownames(p.het) <- c("Liberals","Conservatives","Muslims","People on welfare","Atheists","Feminists")
colnames(p.het) <- c("No logo vs. UC Berkeley","No logo vs. Liberty U","UC Berkeley vs. Liberty U","Joint")


p.het[1,1:3] <- 

c(
if(t.stat.1.1v2.h > 0) round(mean(t.stat.rand.het[,1] > t.stat.1.1v2.h),2) else "*",
if(t.stat.1.1v3.h > 0) round(mean(t.stat.rand.het[,2] > t.stat.1.1v3.h),2) else "*",
if(t.stat.1.2v3.h > 0) round(mean(t.stat.rand.het[,3] > t.stat.1.2v3.h),2) else "*")

p.het[2,1:3] <- 

c(
if(t.stat.2.1v2.h > 0) round(mean(t.stat.rand.het[,4] > t.stat.2.1v2.h),2) else "*",
if(t.stat.2.1v3.h > 0) round(mean(t.stat.rand.het[,5] > t.stat.2.1v3.h),2) else "*",
if(t.stat.2.2v3.h > 0) round(mean(t.stat.rand.het[,6] > t.stat.2.2v3.h),2) else "*")

p.het[3,1:3] <- 

c(
if(t.stat.3.1v2.h > 0) round(mean(t.stat.rand.het[,7] > t.stat.3.1v2.h),2) else "*",
if(t.stat.3.1v3.h > 0) round(mean(t.stat.rand.het[,8] > t.stat.3.1v3.h),2) else "*",
if(t.stat.3.2v3.h > 0) round(mean(t.stat.rand.het[,9] > t.stat.3.2v3.h),2) else "*")

p.het[4,1:3] <- 

c(
if(t.stat.4.1v2.h > 0) round(mean(t.stat.rand.het[,10] > t.stat.4.1v2.h),2) else "*", 
if(t.stat.4.1v3.h > 0) round(mean(t.stat.rand.het[,11] > t.stat.4.1v3.h),2) else "*",
if(t.stat.4.2v3.h > 0) round(mean(t.stat.rand.het[,12] > t.stat.4.2v3.h),2) else "*")

p.het[5,1:3] <- 

c(
if(t.stat.5.1v2.h > 0) round(mean(t.stat.rand.het[,13] > t.stat.5.1v2.h),2) else "*",
if(t.stat.5.1v3.h > 0) round(mean(t.stat.rand.het[,14] > t.stat.5.1v3.h),2) else "*",
if(t.stat.5.2v3.h > 0) round(mean(t.stat.rand.het[,15] > t.stat.5.2v3.h),2) else "*")

p.het[6,1:3] <- 

c(
if(t.stat.6.1v2.h > 0) round(mean(t.stat.rand.het[,16] > t.stat.6.1v2.h),2) else "*",
if(t.stat.6.1v3.h > 0) round(mean(t.stat.rand.het[,17] > t.stat.6.1v3.h),2) else "*",
if(t.stat.6.2v3.h > 0) round(mean(t.stat.rand.het[,18] > t.stat.6.2v3.h),2) else "*")


##  p-values for family-wise comparisons (within outcomes)
p.het[,4] <- 

round(
c(
mean(apply(t.stat.rand.het[,1:3],1,max) 	> max(c(t.stat.1.1v2.h,t.stat.1.1v3.h,t.stat.1.2v3.h))),
mean(apply(t.stat.rand.het[,4:6],1,max) 	> max(c(t.stat.2.1v2.h,t.stat.2.1v3.h,t.stat.2.2v3.h))),
mean(apply(t.stat.rand.het[,7:9],1,max) 	> max(c(t.stat.3.1v2.h,t.stat.3.1v3.h,t.stat.3.2v3.h))),
mean(apply(t.stat.rand.het[,10:12],1,max) 	> max(c(t.stat.4.1v2.h,t.stat.4.1v3.h,t.stat.4.2v3.h))),
mean(apply(t.stat.rand.het[,13:15],1,max) 	> max(c(t.stat.5.1v2.h,t.stat.5.1v3.h,t.stat.5.2v3.h))),
mean(apply(t.stat.rand.het[,16:18],1,max) 	> max(c(t.stat.6.1v2.h,t.stat.6.1v3.h,t.stat.6.2v3.h)))),2)

p.het[2,4] <- "*" # all test statistics have the wrong sign for conservatives


library(Hmisc)
latex(p.het, file = "")



##
##
##  Additional analysis of treatment effect heterogeneity prompted by R3.
##  We treat the liberal and conservative scales are proper covariates and
##  look for treatment effect heterogeneity for the other 4 temperature scales.
library(lmtest)
library(sandwich)

res1 <- matrix(NA,24,3)
dat$affective <- dat$liberals-dat$conservatives


##  Muslims
lm.out.muslim <- lm(muslims ~ as.factor(Tr)*affective, data = dat)
summary(lm.out.muslim)

res1[1:6,] <- coeftest(lm.out.muslim, vcov = vcovHC(lm.out.muslim))[,c(1,2,4)]
rm(lm.out.muslim)


##  People on Welfare
lm.out.welfare <- lm(welfare ~ as.factor(Tr)*affective, data = dat)
summary(lm.out.welfare)

res1[7:12,] <- coeftest(lm.out.welfare, vcov = vcovHC(lm.out.welfare))[,c(1,2,4)]
rm(lm.out.welfare)


##  Atheists
lm.out.atheists <- lm(atheists ~ as.factor(Tr)*affective, data = dat)
summary(lm.out.atheists)

res1[13:18,] <- coeftest(lm.out.atheists, vcov = vcovHC(lm.out.atheists))[,c(1,2,4)]
rm(lm.out.atheists)


##  Feminists
lm.out.feminists <- lm(feminists ~ as.factor(Tr)*affective, data = dat)
summary(lm.out.feminists)

res1[19:24,] <- coeftest(lm.out.feminists, vcov = vcovHC(lm.out.feminists))[,c(1,2,4)]
rm(lm.out.feminists)


res1 <- round(res1,2)
colnames(res1) <- c("coef","SE","p")
rownames(res1) <- rep(c("Intercept", "UC Berkeley", "Liberty U", "Polarization",
					"UC Berkeley x Polarization", "Liberty U x Polarization"),4)

library(Hmisc)
latex(res1, file = "")













.
